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Abstract. In a system of ac-driven condensed bosons we study a new type of 
Josephson effect occurring between states sharing the same region of space and the 
same internal atom structure. We first develop a technique to calculate the long time 
dynamics of a driven interacting many-body system. For resonant frequencies, this 
dynamics can be shown to derive from an effective time-independent Hamiltonian 
which is expressed in terms of standard creation and annihilation operators. Within 
the subspace of resonant states, and if the undriven states are plane waves, a locally 
repulsive interaction between bosons translates into an effective attraction. We apply 
the method to study the effect of interactions on the coherent ratchet current of an 
asymmetrically driven boson system. We find a wealth of dynamical regimes which 
includes Rabi oscillations, self-trapping, and chaotic behavior. In the latter case, a 
full many-body calculation deviates from the mean-field results by predicting large 
quantum fluctuations of the relative particle number. 
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1. Introduction 

The Josephson effect (JE) is a fundamental quantum phenomenon which refiects the 
coherent occupation of two or a few single-particle states by a macroscopic number of 
bosons. It was first predicted [1] and observed [2] for superconductors, and seen later 
in superfiuids [3]. In Bose-Einstein condensates (BECs) one may refer to the internal 
or external JE, depending on whether the two coherently connected states occupy the 
same region of space with different internal atomic states, or two different regions with 
the same spin state. Both the internal [4, 5] and the external [6] JE have been observed. 
The internal JE holds the promise of generating highly entangled quantum many-body 
systems [7]. Here we report on a new type of JE which involves the macroscopic coherent 
occupation of Floquet states in ac-driven BECs. This orbital Josephson ejfect is properly 
neither external nor internal, since the connected Floquet states occupy the same region 
of space with the same internal atom structure, their sole difference residing in their 
orbital state. 

The Josephson effect requires interactions in order to display truly collective 
behavior associated with the macroscopic occupation of two states [8, 9]. In the absence 
of interactions, the resulting Rabi dynamics is merely an amplification of the dynamics 
undergone by a single atom [4, 9]. The role of interactions in ac-driven many-body 
systems is difficult to treat, and frequently simplifications are made such as the two-mode 
approximation [10], mean field theory [11, 12, 13], and its first [14, 15] and second [16, 17] 
order corrections, or the use of an effective description [18] in terms of a static many- 
body Hamiltonian with renormalized parameters. This latter approach is valid when 
the driving frequency is the dominant energy scale. Alternatively such systems can 
be numerically studied by the exact simulation of small clusters [18], or by recently 
developed techniques such as the time-dependent density-matrix renormalization group 
(t-DMRG) method [19], or multiconfigurational time-dependent Hartree for bosons 
(MCTDHB) [20], both of which are able to treat larger systems. Here we develop a 
description of the coarse-grained dynamics of the quantum field operator in an ac-driven 
many-boson system, thus going beyond a mean-field treatment and its first corrections. 
We find that, in the case of resonantly connected Floquet states, the long-time evolution 
can be described by conventional quantum dynamics, by which we mean one deriving 
from an effective interacting Hamiltonian. If the driving provides the only external 
potential, we find that a repulsive interaction in real space translates into an attractive 
interaction in the truncated Hilbert space of a discrete number of resonantly connected 
Floquet states. We apply these findings to understand in depth the role of interactions 
in asymmetrically driven BECs which exhibit the quantum ratchet effect [11, 12, 21, 22], 
together with a wealth of dynamic regimes. 
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Figure 1. Schematics of the orbital Josephson effect, (a) Resonant driving induces 
coupling between the zero-momentum and two opposite momentum eigenmodes. The 
resulting Josephson link between these orbitals, with Rabi frequencies 2r-|-, has a 
negative effective interaction energy, (b) A possible departure from exact resonance 
lifts the degeneracy of the quasi-encrgics. 



2. (t,t') formalism in second quantization 

Our starting point is the equation of motion for the field operator of a system of 
interacting bosons, 

idtipix, t) = [H{x, t) + \'iij^]'iij{x, t), (1) 

where H{x,t) = —\dxx + V{x,t) includes a time-dependent potential, and A is the 
effective strength of the contact interaction in a quasi-one-dimensional ring of radius R. 
We set ^ = 1 and measure all energies and frequencies in units of h'^/MR'^, with M the 
atomic mass. Furthermore, our unit of length is R, so that the circumference of the 
ring is 27r. Such rings, as schematically depicted in Fig. la, have been proposed and 
experimentally realized - see e.g. Refs. [23, 24], among many other ring experiments. 
A possible particle flow can be observed via standard time-of-fiight measurement 
techniques [23]. The field operator i^{x,t) obeys periodic boundary conditions in x, 
and its algebraic structure is given by the standard equal-time bosonic commutation 
relations: [ip{xi,t),ip'^ {x2,t)] = 6{xi — X2) and [ip{xi,t),ip{x2,t)] = 0. 

In order to study the long-term dynamics of the many-body system, we generalize 
the {t, t') formalism to treat second-quantized operators. This formalism was originally 
developed for single particles [25], and later adapted to treat helium [26] always within 
first quantization. We do this extension by introducing an additional parameter t' 
which the field operators depend on, and impose periodic boundary conditions in t' 
with period T 

i){x,t' + T-t)=i){x,t']t) (2) 
We define ^/'(a;, t'] t) as the solution of the equation of motion 

i9fV'(x, t'] t) = [H{x, t') - idf + X'ip^]'ip{x, t'] t) , (3) 
together with the initial condition 

^(x,t';0) = ^(x,0). (4) 

The periodic boundary conditions in t\ which are imposed initially via Eq. (4), are 
guaranteed to be preserved over time t due to the periodicity of H{x,t'). 

We note that the field operator ip{x,t']t) acts on the same Fock space as the 
Heisenberg field operators ■ip{x,t), for all t' and all times t. 
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We may choose t' to be restricted to be a function of t given by t'^{t) = t + t, 
with r G [0,T]. Then the resulting operator i^T-{x,t) = ijj{x,t'^;t) = ■ip{x,t + r; t) is 
the solution of the Heisenberg equation (1), but with a shifted switching HT-{x,t) = 
H{x,t + r) and initial condition ipr{x,0) = ip{x,0). This can be shown by calculating 
the total time derivative along the line t' = t'^{t): 



idtiprix, t) = ^^^(a;, K; t) 



. dip dt'^ , dip 

^dElF^^~dt 



t'=t+T 



t'^=t+T 

= [H{x,t + T) + \iP^]tP, (5) 

where we have used Eq. (3) and the identity dt'^/dt = 1. In other words, ip{x,t'^;t) = 
iprix, t) is but a solution of the Heisenberg equation of motion for a Hamiltonian shifted 
in time by r with the initial condition ipr{x,0) = ip{x,0). In particular t' = t (r = 0) 
gives the solution to Eq. (1). This implies the following commutation relations for the 
solutions of Eqs. (3)-(4) 

['ip{xi,t';t),i>\x2,t';t)] = 6{xi- X2), 

(6) 

[tP{xi,t';t),tP{x2,t';t)] = . 

The advantage of the {t, t') approach is that it provides a natural separation of 
timescales. The t' coordinate describes the behavior of the system on timescales shorter 
than the driving period T, while the long-time dynamics is described by t - see Eqs. (7)- 
(8) below. Furthermore, the originally time-dependent problem [see Eq. (1)] is mapped 
to a formally time-independent one [see Eq. (3)], which opens the possibility of using 
resolution methods developed for such type of problems [21, 25, 26]. 

The field-operator extension of the (t, t')-formalism which we have presented is 
rather general. It is neither restricted to bosonic systems, nor to the ring geometry we 
study here. It can be straightforwardly translated to the case of fermionic or mixed 
systems, as well as to other trap-geometries. 

Due to the periodic boundary conditions in t', it is natural to express the t'- 
dependence of xjj{x,t';t) via a Fourier decomposition, and we may approximate the 
dynamics by incorporating only a few modes. 

To include the effect of a weak time-periodic perturbation, V{x,t) = V{x,t + T), 
it is convenient to work in the representation of unperturbed {V = 0) Floquet states. 
The first-quantized version of (3) was studied in Ref. [21] in the absence of interactions. 
The unperturbed stationary states evolve as 

ipem{x, t'; t) = (pem{x, t') exp{-ie^g^t) , (7) 

where = — mu and £, m are integers labeling the Fourier modes 

4>em{x, t') = exp{iix — imut') , (8) 



/27r 

which are the unperturbed Floquet states with quasi-energy e^^, with u = 2ii/T. The 
operators in this representation, 

a^m{t) = 1; [fdx dt'(P;^{x, t'; t) , (9) 



T 
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satisfy the equation of motion 

(10) 

Zn ^ — ^ 

e'm'i"m" 

and the commutation relations (6) translate to 
''^^[0'e,m'+mit),al,^,(t)] = See'Smo , 

(11) 

^[ai^m'-m{t),aem'{t)] = , for all m. 

m' 

If one were to view m as just an additional standard (orbital or spin) quantum number, 
these commutation relations would appear unconventional (very much as Eqs. (6) would 
look unconventional if t' were regarded as an extra space variable). However, Eqs. (11) 
follow naturally from the Fourier transformation of Eqs. (6). 

Here Vim^i'm' is the matrix element of the driving operator between two unperturbed 
states, as given in Eq. (8). When the driving frequency is such that the system is at 
or close to resonance, only a few states (all with the same or similar value of e^^) are 
relevant [11, 12, 21]. Importantly, in that subspace the index m is uniquely determined 
by £. Thus within that truncated space we can drop the index m. As a result, Eqs. (11) 
become equivalent to the standard bosonic commutation relations [ai{t),dl,{t)] = 5u' 
and [a lit), a I' it)] = 0. 

3. Resonant driving 

If we calculate the effective matrix elements connecting the various resonant unperturbed 
Floquet states, we are left with a conventional few-mode boson problem whose dynamics, 
given by Eq. (10), can be studied using established techniques. The resulting dynamics 
between Floquet states reflects the coarse-grained, long-time dynamics of the true 
quantum state evolution. 

Within the degenerate (or almost degenerate) subspace where conventional 
commutation relations apply, Eq. (10) can be viewed as a Heisenberg equation deriving 
from the interacting many-body Hamiltonian 

H = ^2 ^l^f- + X/ r«'aja£/ + ^ ^ ^ a\^f^„a\,_^„aiaii , (12) 



47r ^ ^ 

«' a' I" 



where ei = mod u, Tui may allow for second-order processes between resonant states 
mediated by a non-resonant state [21], hi = a\ai is the occupation of state £, and can 
only take values or f-£ ii I ^ f , and if £ = f . 

We remark that the derivation of an effective Hamiltonian does not require £ to 
describe a plane wave. In such a case, the interaction term in Eq. (12) would contain a 
sum over four orbital indices, and the matrix elements would depend on all four involved 
orbitals. 
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The conservation of total particle number he = N, together with the identity 
^n,n,, = iV2-^n2, (13) 

leads to the interesting result that, up to a constant (2A^^ — A^)A/47r, the 
Hamiltonian (12) is equivalent to 

E = ^e\he^^Yu'a\ae - ^^h\ . (14) 

t w I 

Thus a repulsive interaction in real space translates into an attractive interaction in 
(angular) momentum space. The resulting energy gain in the macroscopic occupation of 
a single state in gases with repulsive interactions underlies the stability of Bose-Einstein 
condensation, as discussed in Ref. [27]. 

To derive (14), we have used the selection rules (stemming from the simultaneous 
conservation of momentum and quasi-energy) which are implicit in Eqs. (10)-(12). 
Interestingly, in those dynamical regimes where the truncated resonant space is 
effectively reduced to two modes, Eq. (14) applies for orbitals of arbitrary shape (not 
necessarily plane waves), but with an interaction strength that is state dependent. 

When a few modes are macroscopically occupied, Eq. (12) describes a Josephson- 
type link between states ii' with Rabi frequency 2ra/. Remarkably, this Josephson link 
takes place between atom states with the same internal state and occupying the same 
region of space. Thus it is not appropriate to describe it as an internal or external JE. 
The Hamiltonian (14) describes what can be termed the orbital Josephson effect, which 
here occurs between resonant Floquet states in ac-driven Bose condensates. 



4. Application to the quantum ratchet 

As a particular application of the method developed here to treat the many-body 
problem in driven systems, we consider a BEC subject to the asymmetric driving 

V{x, t) = K[sm{kx) + a sm{2kx + (f)] [sin(wt) + /3 sin(2wt)] , (15) 

which has been numerically studied [12] and experimentally implemented by Salger et 
al. [28] on an extended optical lattice. Since the driving in Ref. [28] conserves quasi 
momentum, the here considered periodic boundary conditions can be related to the 
experimentally more convenient lattice system. For ip ^ ^ir, |7r it yields a coherent 
quantum ratchet provided that both a and /3 are nonzero. 

We consider a system initially prepared in |£) = |0). It was analytically shown in 
Ref. [21] that, if we drive it with resonant frequency uj = 1 and small amplitude K, 
the subsequent evolution will mix the initial state only with the states \i) = | ± 2). 
These three states all satisfy the resonance condition = 0, and we label them with 
the indices ±,0. The resulting three- level system Hamiltonian is 

-^3LS = r+d^_ao + r„aLao + H.c. - ^ hi - 2A(n+ + n_), (16) 



Orbital Josephson effect and interactions in driven atom condensates on a ring 



7 



where i' takes values ±, 0, and A = u — 1 accounts for a possible small detuning that 
shifts the quasi-energies in Eq. (14) to = — (£^/2)A. The resulting level structure is 
schematically depicted in Fig. lb. Near resonance the tunneling parameters T± can be 
calculated analytically [21]. For A ^ 1, one finds 



This result demonstrates that the quantum ratchet current (which requires |r+| ^ |r_|) 
originates in the interference between first- and second-order processes in the driving 
strength [21]. In Eq. (17) we have neglected the possible effect of interactions on the 
second-order tunneling terms. Without interactions, the initial state |0) couples only 
to a state \a) which is an asymmetric combination of |±), yielding an average current 
which is half that carried by \a) [21]. It has been noticed that interactions destroy the 
coherent current [12]. Below we study the role of interactions in depth and show that 
they yield a rich variety of dynamical regimes. 

5. Numerical results 

We employ three different approximations to investigate the many-body problem: full 
numerical resolution of the time-dependent Gross-Pitaevskii equation (FGP) using a 
fourth-order Runge-Kutta method with a fixed time-step; study of the GP equation 
in the truncated space of the three resonant states 0, ± (3GP); and resolution of the 
many-body problem in the three-level system (3LS), via exact diagonalization of the 
effective interacting Hamiltonian (16) with up to 40 atoms. Importantly we note that 
3GP can be obtained both as a truncation of FGP or a mean-field version of 3LS. We 
always assume the condensate to be initially in state 0. 

5.1. Time- dependent current 

Figure 2 shows the expectation value of the ratchet current {I(t)) - given by the mean 
momentum per particle - over the first 200 driving cycles. The regimes of weak, 
moderate, and strong interaction [as characterized hj g = X{N — 1)], are all calculated 
with the three methods described above. For weak interactions, the two truncated- 
space calculations (3GP and 3LS) yield similar Rabi oscillations, both differing from 
FGP in that they do not display high-frequency dynamics, as revealed both in the 
Fourier spectrum and in the time dependence of the current. The fast dynamics of FGP 
disappears if the stroboscopic current is plotted, as shown in the shaded region. 

For strong interactions, the system tends to remain trapped in the initial state |, 
as may be expected since Eq. (16) has the structure of a Josephson link with negative 
interaction energy [29]. The three numerical methods predict differing behaviors of the 
resulting small current, but agree on the position and strength of the main peak in the 




(17) 



I In the absence of decoherence the initial state would mix with resonant states ± after a very long 
time. 
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Figure 2. Time evolution and power spectrum of the particle current. Three 
regimes are depicted: (a,b) weak interaction strength (g = 0.01), showing regular Rabi 
oscillations; (c,d) intermediate interaction ((7 = 0.08), displaying chaotic dynamics; and 
(e,f) strong interaction (g = 0.5), showing self-trapping. Parameters: K = 0.2, ip ^ 0, 
ui = 1 (A = 0),a=:/3 = 0.2. In the shaded region we present the stroboscopic 
current, which erases the short time dynamics responsible for the peaks in the white 
region of the current spectrum. Three different calculation methods are employed: Full 
Gross-Pitaevskii (FGP), three- level mean- field description (3GP), three- level many- 
body description (3LS) (see section 5). In (a,b) the green curve almost coincides with 
the red curve. 



spectrum. For intermediate interaction, discrepancies between the three methods soon 
appear, which suggests chaotic behavior. This is confirmed and studied further below. 

Within the accuracy hmited by the w-samphng, the positions of the main peaks 
coincide exactly in the Rabi regime, and vary by about 10% in the self-trapping case. 
The peak heights vary by about 6% in the Rabi case, and about 9% between FGP and 
3GP in the the self-trapped regime. However for 3LS, the discrepancy with FGP (and 
3GP) is worse (differing by a factor of 3), which coincides with the observation that for 
the self-trapped dynamics one has to consider particle numbers far beyond = 40 in 
order to get an agreement between the full many-body description and the mean-field 
result [30]. 
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Figure 3. (a) Continuous time average of the quantum ratchet current at long times 
(400 driving cycles) for different amplitudes and interaction strength, (b) Third largest 
eigenvalue of the time-averaged one-particle density matrix after the same interval. 
Parameters: w = l(A = 0),a = /3 = 0.2, = 0. 

5.2. Time- averaged current 

Figure 3a shows the continuous time average of the ratchet current 

m - \ 

after 400 cycles, computed within 3GP, as a function of g and K at exact resonance. 
We assume |r_|_| > |r_|; see Eq. (17). For small i^, the perturbative scheme of [21] 
applies, and at weak interaction a nonzero ratchet current results from Rabi oscillations 
between |0) and \a). For stronger interactions the ratchet current disappears, because 
the system remains in the initial state |0) due to the large negative energy associated with 
the initial macroscopic occupation of |0) [see Eq. (14)]. This can be viewed as a form 
of macroscopic quantum self-trapping [31], albeit with a negative effective interaction 
energy. We find that just below the critical interaction strength the system oscillates 
between |0) and |+). This is consistent with the fact that the threshold interaction value 
obtained from that assumption, g^ = 87r|r+|, agrees well with the numerical simulation 
(see solid curve in Fig 3a). Finally, we notice the existence of a region of suppressed or 
weakly reversed current for stronger driving and intermediate interaction. 



dt' {I{t')) 



5. 3. Effect of detuning 

Figure 4a shows the effect of departing from exact resonance. In general, the effect 
of A > is that of shifting the behavior of current towards higher interactions. In 
particular we note that, starting from the white-blue region of suppressed or reversed 
current, interactions tend to restore the positive ratchet current. This effect can be 
understood if one notes that, for positive detuning, the effective degeneracy between 
the three states ±,0 disappears because |0) acquires a higher energy. Degeneracy is 
restored due to the massive initial occupation of |0), which lowers its energy via the 
attractive mean-field interaction in Eq. (14). It thus becomes clear that the effect of a 
positive detuning is counteracted by an increase of the effective attraction. In particular. 
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Figure 4. (a) Quantum ratchet current as a function of the interaction strength and 
the detuning. The behavior along the dashed line {g = 0.1) is shown in (b) for four 
different calculations (3GP is computed at two different times, 400T and lO^T). The 
black dots in (b) show the maximal Lyapunov exponents and refer to the scale on the 
right (rip; is the Rabi frequency for the non- interacting case, see Ref. [21]). Poincare 
sections are shown for the representative points (c-f), obtained from points in the 
hyperplane =0.2. The plotted variables are the instantaneous ratchet current 

(/) and the relative phase (p+ between states |0) and |+). Parameters: a = /3 = 0.2, 
K = 0.2, ip = 0. In (a-b) averages are performed over 400 cycles unless otherwise 
indicated. 

the starting degeneracy recovered with the help of interactions permits the onset of Rabi 
oscillations known to be essential for the emergence of a ratchet current [21]. 

5.4- Phase dependence and interaction symmetry 

In Fig. 5 we show how the the ratchet current depends on the choice of initial driving 
phase. A change </) = — tt corresponds to K ^ —K and thus to a reversal of the 
current. On the other hand, we note that for ip = ^n, |7r the ratchet current is always 
zero, because parity is restored for those values of ip. Figure 5 also shows that the 
long-time averaged current (but not the instantaneous one) is an even function of the 
interaction strength. 
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Figure 5. Current dependence on ip. Ratchet current in 3GP as a function of the 
phase (p and the interaction strength g. Parameters: a; = 1 {A — 0), K = 0.2, 
a — f3 ^ 0.2. Averages are performed over 400 cycles. 



5.5. Chaotic dynamics 

Figure 3b shows the population of the third eigenstate of the time-averaged one- 
particle density matrix (calculated in 3GP) 



where A^{t) (/i = 0, ±) are the time- dependent expansion coefficients of the condensate 
orbital. Comparison with Fig. 3a shows a clear correlation between the occupation 
of a third state and the suppression or weak reversal of the ratchet current. The 
dynamics of a macroscopic condensate in a three-level system is that of two coupled 
non-rigid pendula, which beyond the harmonic regime can be expected to be chaotic. 
This behavior is confirmed from the inspection of Figs. 4b-f. In Fig. 4b we show a plot 
of / as a function of A for g = 0.1, obtained with the three different calculation methods 
described before. For the values of A marked with labels c-f, we show in Figs. 4c-f the 
corresponding Poincare cross sections. The one point which shows chaotic behavior is 
e, which falls in the white-blue region of Fig. 4a. We have just noted (see Fig. 3a) that 
such a region is correlated with the occupation of a third eigenstate and thus with the 
possible emergence of chaotic behavior on the coarse-grained time scale. This underlines 
the connection between three-level occupation, chaotic behavior, and current reduction 
or reversal. A calculation of the maximal Lyapunov exponent Am (given by black dots in 
Fig. 4b) further confirms the chaotic character of this parameter region, as Am acquires 
positive finite values there. For more details, see Appendix A. 

Figure 4b also indicates that in the chaotic regime the convergence in time is slower 
than in the regular regime, as can be seen by the discrepancy between 3GP calculations 
performed over 400 or 10^ driving periods. We emphasize that the current reversal 
feature in the chaotic region is not just a transient but refiects a lasting behavior; we 
have checked it up to 10^ periods in some cases. Finally, we note that all calculations 
performed in the FGP approximation indicate that, for large K and g close to the 
self-trapping transition g^ (beyond which the system essentially remains in |0)), the 
behavior undergoes a significant infiuence from states beyond the three-state basis 0, ± 
(not shown). 




(19) 
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Figure 6. Correlated occupation in the chaotic regime, (a) Normalized occupation 
numbers = [uij/N and (b-c) normalized particle number correlations Cij of 
the (instantaneous) natural orbitals (cigcnstatcs of the reduced one-particle density 
matrix) for iV = 40 atoms. The time-point imax of the local maximum of Coo [marked 
as a blue square in (c)] for different particle numbers TV is given in (d) and its value 
CJ^^ is shown in (e) - same for the maximum of Cn (red triangles) and the minimum 
of Coi (black circles). In (a-c) parameters are the same as in Fig. 2c-d. In (d-e), the 
parameters are also those of Fig. 2c-d except for the total particle number N ^ which 
is here plotted as a variable. 



5. 6. Two-particle correlation functions 

It was noted in Ref. [32] and hinted in a previous work [14] that chaotic motion within 
GP indicates the end of vahdity of this approximation. This was recently reexamined in 
Refs. [33, 34]. Within the 3LS, we are able to test the validity of GP directly for relatively 
high particle numbers. We do this by analyzing the relative occupation numbers of the 
natural orbitals, which are those obtained by diagonalizing the reduced one-particle 
density matrix, defined in a many-body context as 

= (20) 

For the two limiting cases of Rabi oscillations and self-trapping, the relative occupation 
of the condensate orbital (defined here as the most occupied eigenstate of p^^^) stays 
above 98% during the first 400 driving cycles, which justifies a mean-field treatment. In 
contrast, for the chaotic regime, we find that the occupation of the condensate orbital 
decreases abruptly after some 120 cycles (about 2 Rabi periods) and a second natural 
orbital gets macroscopically occupied - see Fig. 6a. We find a low sensitivity of this 
behavior to the particle number. 
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In this chaotic regime, the macroscopic occupation of more than one orbital does 
not occur via a simple fragmentation of the condensate. This can be seen in the particle 
number correlation functions Cij = (niUj) — {ni){nj), where i,j refer to the natural 
orbitals. In a fragmented condensate, Cij = ACij/N^ would be zero (for oo), while 

for a noon (cat) type state such as e.g. 

|noon) = ^(|A^00) + |0A^0)) (21) 

its absolute value approaches unity, reflecting the creation of macroscopic particle 
number correlations. Figures 6b-c show the time-evolution of Cij for the same parameter 
set as in Figs. 2c-d. We see that, with the decrease of the condensate fraction, some 
particle number correlations increase in a non-negligible way. In order to check the 
scaling of this behavior with the particle number, we focus on the local maximum 
of the variance of the condensate occupation (Coo) for different particle numbers. 
Figure 6d shows that this maximum gets shifted to later times for higher A^. The 
value of this maximum stays above 0.2 for all considered particle numbers and increases 
monotonically for > 25, as shown in Fig. 6e. This indicates a scaling of the 
correlations Cij with A^^, which in turn reveals noon-like (21) behavior, i.e., a dynamics 
dominated by a few many-body configurations differing by macroscopically large relative 
particle numbers. Similar observations have been made by Weiss and Teichmann [35], 
investigating a driven double-well system. 

Thus we find that in this chaotic regime more than one orbital gets macroscopically 
occupied. At the same time the particle number correlation functions scale with the 
square of the total particle number. This indicates a complex dynamics in which the 
macroscopically large particle number correlations become a substantial feature. 

6. Conclusions 

In summary, we have introduced a formalism to account for interactions in ac-driven 
many-body systems. In the case where only resonant states effectively intervene, we 
derive a Hamiltonian involving conventional operators which describes the long time 
dynamics of the driven system. Although we have focused on the case where the undriven 
states are plane waves, the method we have developed applies to an arbitrary set of states 
that become resonant in the presence of ac-driving. We have applied the method to the 
calculation of the coherent ratchet current carried by an asymmetrically driven atomic 
condensate, and compared it with continuum and truncated mean-field descriptions. 
We have found a rich dynamical behavior with crossovers from self-trapping to chaotic 
behavior to regular oscillations. In the latter case, we find a strong departure from 
the mean-field picture and in particular an appreciable increase of particle number 
correlations. 
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Appendix A. Calculation of Lyapunov exponents 

In this appendix we present details on the calculation of the maximal Lyapunov exponent 
Am given in Fig. 4b. The Lyapunov exponent reflects how fast the distance between two 
neighboring trajectories grows over time. We quantify the distance between two Hilbert 
space vectors via the Euclidean norm: 

d2 i^l ,i^2) = -^2|^1 -^2), (A. 1) 

with — 4^2) = I'ipi) — 14^2)- The maximal Lyapunov exponent for the trajectory |^/'o(^)) 
is defined over the double limit 



Am = lim lim - log 



d{ij^{t),Mt)) 



(A2) 



.rf(^i(0),^o(0)). 

where |'i^o(^)) and are both solutions of the CP equation. Because of the 

Hamiltonian structure of the nonlinear Schrodinger equation, it is convenient to allow 
only perturbations that yield states with the same value of the Hamiltonian [33]. Note 
that in our effective description the Hamiltonian reflects the system's quasi-energy 
instead of its energy. Within 3GP, the system's dynamics is given by the 3-tuple 
A{t) = Ao(t), reflecting the expansion coefficients of the condensate 

orbital with respect to the three modes ±, 0. In order to extract the Lyapunov exponent 
from the system's dynamics, we consider the time evolution of the two initial states 

As+ = Ns{i5,l,0) and A^^ = N5(0, 1, (A.3) 

where N^ = 1/Vl + 5^ is the normalization factor. The parameter 5 parameterizes the 
initial distance, \/2S, between the trajectories As±{t). The Hamiltonian of both the 
states A5-1- has the same value, and for S = 0, we recover the initial state considered by 
us. We denote their time-dependent separation by ds{t) = d2{As+{t), As-{t)), and 
use this quantity to estimate the Lyapunov exponent. The limit tpi ipQ in the 
definition (A. 2) translates into the limit 5 — 0. 

Figure Al shows the time evolution of ds{t) on a logarithmic scale for chaotic (a), 
and regular dynamics (b). Two different values of 5 are considered. For chaotic dynamics 
the distance between the two trajectories ds{t) grows over several orders of magnitude, 
until it saturates at a value around unity. This saturation occurs because the distance 
ds{t) is bounded from above by the value 2, since both states As± are normalized to 
unity. The value of Am is estimated by fitting an exponential to the rising slope of 
dsit) before saturation is reached. For the case shown we obtain Am ~ 0.61 (0.51) for 
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Figure Al. Trajectory distances over time. Two paradigmatic cases are considered: 
(a) chaotic, and (b) regular dynamics. After 3000 T (vertical line), the time 
propagation is inverted to check the accuracy. In (a), a linear (on a log-scale) fit 
is included for both values of S, given by the straight gray lines. The parameters are 
the same as in Figs. 4c (b) and 4e (a). 



5 = 10~^ (IQ-^). The estimated values thus have a relative difference of about 18%, 
which implies a high relative uncertainty. Nevertheless, the obtained value for Am is 
clearly positive, indicating chaotic dynamics. 

In the regular regime ds(t) deviates by less than 10% of its initial value during the 
depicted time range, for both values of S. Chaos can be precluded when the global 
maximum of ds{t) can be decreased by choosing smaller values for S. 

In all cases, the accuracy of the simulation was checked by propagating both the 
states A^-i- backward in time after t = 3000 T. The simulation is reliable, when ds{t) 
comes back to its starting value after t = 6000 T. This check is necessary to ensure that 
the observed separation of trajectories is not due merely to numerical errors. 

We additionally calculated the maximal Lyapunov exponent by a different method 
described e.g. in Ref. [36]. The results obtained by that approach confirm the results 
presented here, by identifying the same parameter regions as chaotic (and regular) and 
yielding values Am of similar magnitude. 
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